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VIBRATION AND BUCKLING OF PRESTRESSED SHELLS OF REVOLUTION 1 

By Paul A. Cooper 
Langley Research Center 

SUMMARY 

A linearized set of equations is developed for the infinitesimal vibration and 
buckling of an axi symmetric ally prestressed thin shell with an arbitrary meridional con- 
figuration. A finite -difference numerical procedure is given for finding the natural fre- 
quencies, buckling loads, and associated mode shapes, and the procedure is applied to 
calculation of natural frequencies of an unstressed, simply supported cylinder. A 
closed-form solution is obtained for the simply supported cylinder and is used to verify 
the numerical procedure, which is then used to solve some other example problems 
including the vibrations of particular shells of zero and positive Gaussian curvature. 

INTRODUCTION 

In the design of shell structures for launch vehicles, planetary atmospheric entry 
probes, or similar structures, knowledge of the natural frequencies and mode shapes of 
the systems is of fundamental importance in determining their dynamic behavior. For 
strength considerations of such shells, knowledge of the structural stability of the con- 
figurations under varying aerodynamic pressure distributions is important. These shell 
configurations are often too complex to allow solutions in closed form, and numerical 
techniques are appropriate for investigating their dynamic and static behavior. 

Several numerical methods have recently been developed for static stress analysis 
and non-prestressed free vibration analysis of general shells of revolution. Static stress 
analysis is considered in references 1 to 4 and non-prestressed vibrations are considered 
in references 5 and 6. Membrane and flexural vibrations of toroidal shells are treated 
in references 7 and 8 on the basis of the numerical approach of reference 1. 

The investigation of prestressed vibrations and buckling of shells requires a con- 
sideration of the nonlinear shell equations. The object of this paper is to develop 

1 The information presented herein is to be included with additional material to be 
offered in partial fulfillment of the thesis requirements for the degree of Doctor of 
Philosophy in Engineering Mechanics, Virginia Polytechnic Institute, Blacksburg, 

Virginia. 
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differential and difference equations that govern the asymmetric vibration and buckling 
of a class of shapes of prestressed shells of revolution. These equations are based on 
the nonlinear theory of reference 9. A finite -difference procedure similar to that of 
reference 1 (utilizing and extending the ideas of refs. 7 and 8) is then formulated to 
obtain solutions of the equations. A comparison with known solutions (refs. 10 to 13) 
for free vibrations and buckling of a simply supported cylinder indicates the validity of 
the equations derived in the present paper. In addition, consideration is given to some 
problems which have not been previously treated in the literature. Natural frequencies 
are calculated for a cylinder with one end simply supported (with in -plane displacements 
free) and the other end clamped. Also, a comparison is made between natural frequen- 
cies for a cylindrical shell and a similar shell with constant positive Gaussian curvature 

SYMBOLS 


a reference length 

B extensional stiffness (eq. (16)) 

D bending stiffness (eq. (16)) 

e x ,e^,E x ,Eg nondime nsional prestress parameters (see eqs. (19) and (52)) 

E Young's modulus of elasticity 

h shell thickness 

i station number, i = 0, 1, 2, . . . N 

k x ,kg nondimensional curvatures (eqs. (12)) 

K critical buckling load parameter 

m number of axial half-waves for cylinder (see eqs. (26)) 

m moment variable (see eqs. (20)) 

modified moment resultants associated with perturbed state 
n number of circumferential waves 
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N 


total number of intervals along the meridian 


NpN^N^N 


N ^ N 0 


r 


R^,Rg 

s 


S 

t 


u,v,w 


q^Q^Qq modified stress resultants associated with perturbed state 

modified prestress stress resultants 
surface loading 

nondimensional radius of cross section, p/a 
principal radii of curvature (see fig. 1) 
total meridional arc length 

total nondimensional meridional arc length, s/a 
time 

displacement variables defining perturbed state (see eqs. (20)) 


U,V,W displacements in meridional (|), circumferential (9), and normal direc- 

tions, respectively, of undeformed middle surface defining perturbed 
state (see eqs. (20)) 


x 


nondimensional meridional coordinate, £/a 


_ 1 dr 
^ r dx 

A 

e £’ e 9’ e £9 

9 


K £> K 9 ,K Zd 


X 

4 

v 


length of interval between stations, S/N 

middle -surface strains associated with perturbed state 

circumferential coordinate in undeformed shell 

middle -surface bending strains associated with the perturbed state 

h 

thickness parameter, — 

a 

Poisson’s ratio 
mass density 


meridional coordinate in undeformed shell 


cross-sectional radius (fig. 1) 


Qlty 




CO 


2 __ a 2 a> 2 y(l - p 2 ) 
E 


sr = 


middle-surface rotations associated with perturbed state 
prestress meridional rotation 
natural frequency 

frequency parameter 


Notations Used to Identify Load and Deformation Variables: 

Unmarked variables indicate variables associated with perturbed state only. 

(~) indicates modified variables associated with the total deformation 

(”) indicates modified variables associated with the prestress state only 

( )° indicates physical stress resultant quantities 

The comma before a subscript denotes differentiation with respect to the following 
subscripted variable. A dot over a symbol indicates differentiation of the quantities with 
respect to time. 

1x4 column matrices: 

Z ^ dependent variable 

4x4 matrices: 


Aj,Bj,Cj,Dq,D 

e jk’ f jk 

F jk’ G jk’ H jk 

P i 

a, 13 


^.,Eq,E^ difference -equation coefficients 

boundary -equation coefficients 
equilibrium -equation coefficients 
recursion coefficients 
boundary-condition-selection matrices 
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DEVELOPMENT OF GOVERNING EQUATIONS 


The shell geometry is illustrated in figure 1. The location of points on the middle 
surface of the shell is described by the principal coordinates (£,0), where £ is the 
meridional distance measured on the 


middle surface from one boundary, 
and 0 is the circumferential angle. 
Since the shell is axisymmetric, it is 
completely described by the merid- 
ional shape parameter p(£) which is 
the radial distance from the axis of 
revolution to the middle surface of the 
shell. 

The principal radii of curvature 
of the middle surface, R^(£) and 
R e (0, are given by: 




Figure 1.- Shell middle-surface geometry. 


The shell is assumed to have a constant thickness h measured along the normal 
to the middle surface and boundaries at £ = 0 and | = s, where s is the total merid- 
ional arc length. The material is assumed homogeneous and isotropic with mass den- 
sity v, Young's modulus of elasticity E, and Poisson's ratio p. 


Governing Nonlinear Equations 

General nonlinear shell equations in which strains are assumed to be small and 
rotations, moderately small, are given in reference 9. For a shell of revolution, these 
equations become, when inertia terms are added, 
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where the comma before a subscript denotes partial differentiation with respect to the 
succeeding subscripted independent variables (£ or 0 ) and dots over a quantity denote 
differentiation with respect to time. The equilibrium equations (2a), (2b), and (2c) rep- 
resent the sum of forces along coordinates of the undeformed surface, and the perturbed 
displacements of the middle surface U, V, and W are measured in the direction of 
the coordinates of the undeformed surface with W measured positive along the outward 
normal. 

The boundary conditions considered on the edges £ = 0 and £ = s may be chosen 
from any combination of the following four pairs of quantities in which either quantity 
(but not both) of each pair is prescribed: 
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N f or U 




«J a w + 3pe + S ») 5 


or V 


^i; + p^£0,0 or 


W 


or 


(3) 


The equations have been derived by use of the Kir chhoff -Love assumptions; that is, 
normals to the undeformed middle surface remain normal to the deformed middle surface, 
normal strain is zero, and the normal stress is negligible. Rotary inertia terms have 
been neglected in the moment equations. Modified stress and moment resultants have 
been used in the development of equations (2) and (3) and are defined as follows: 


N N ! 


Mi 


R, 


n M fl 
N = N 0 - 

"0 "0 r 


0 


M 


N I0 = V 


ZL 


R/ 


M 


N 0£ " N 0£ 




(4) 


= M° 


M e = M° 


' >V = |(V + M 0f°) 


The modified transverse shear stress resultants and Q 0 may be found by applying 

the definitions in equations (4) to the moment equilibrium equations (eqs. (2d) and (2e)). 
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The quantities N. 


N, 


N 




N 


OZ 


Q^°, and 


Q, 


middle -surface stress resultants (fig. 2(a)), and the quantities M^°, Mg°, , 


Normal 



Normal 



(b) Moment resultants and rotations. 
Figure 2.- Middle-surface quantities. 


and M 


represent the total 
M i 0 

represent the 


o 


ez 

total middle -surface moment 
resultants introduced by the 
same combined effect 
(fig. 2(b)). No attempt is made 
to relate these stress and 
moment resultants to the dis- 
tribution of stress through the 
thickness of the shell. The 
equations of reference 9 have 
been derived without depend- 
ence on such a relationship; 
thus, any formulation consist- 
ent with thin-shell theory is 
acceptable. According to ref- 
erence 9, the addition of terms 
like M°/R to the N° quan- 
tities in the stress-strain 
relations does not introduce 
errors any greater than those 
already introduced by neg- 
lecting transverse shear flex- 
ibility in the Kirchh off -Love 
hypothesis. Consequently, the 


N quantities may be treated as stress resultants without introducing an inconsistency in 
the thin-shell analysis. 

The sum of the moments about the normal direction is 

M q ° „ M ta ° 

N 0 - -M- - N ° + — £2- = 0 

W Rg R^ 

Hence, from the definition of the modified stress resultants 


N 0 £ " N |0 


(5) 
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Therefore, the sixth equilibrium equation, that of equilibrium of moments about the nor- 
mal, is identically satisfied by symmetric modified stress resultants. 


Vibration and Buckling Equations 

In the derivation of the vibration and buckling equations, the total rotations and 
total stress and moment resultants are separated into the parts associated with the initial 
axisymmetric prestress and the parts associated with the infinitesimal perturbed dis- 
placements about the prestressed state. The perturbed quantities can be time dependent 
for infinitesimal vibration investigations or static for buckling investigations. Symbols 
with bars represent those quantities associated with the prestress conditions and symbols 
without bars represent those associated with the perturbed state. Thus, the total stress 
state may be completely described by these quantities from the relations: 


W N i 

W N 0 


^0 = N |0 

Q| - Q| + 
^0 = Q 0 


M fc = M t + M t 
? e ? 


M e = M 0 + M e 


M £0 - M |0 


The total rotations are 


<p^ = <?£ + <p% 
V e =<p e 

<p= <p 


(7) 
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and the total surface loading is described by 


P 


I 



P = P 


( 8 ) 


since no additional surface loading is assumed to be associated with the perturbed state. 
Substitution of equations (6), (7), and (8) in equation (2a) yields 


If V ^ i V V 

+ jdf N « + pN l,« + ' d| n 9 + rJ q S + %( r ^ ' 


' rT v ( n ( ' n r ‘ I 


_e_ 


( N ^ + N e) 

«V Ke( s « + S e)} 


+ ) - -2- cp y N 


R, T{ 


(9) 


where 0-derivatives of barred quantities vanish as a result of axisymmetry of the pre- 
stressed state. The prestressed shell is in equilibrium; thus, the sum of the terms 
enclosed by the first set of braces in equation (9) vanishes identically. Furthermore, 
the perturbation of the shell away from the prestressed configuration is governed by 
linear theory. Therefore, the nonlinear terms enclosed by the second pair of braces are 
neglected. 

The term in equation (9) enclosed by the third pair of braces (i.e., the interaction 
between the prestress deformation and the perturbation stress resultants) is usually neg- 
lected in the procedure followed in the classical linearization process for a cylinder. If 
this term is neglected, the general assumption is made that the prestress rotation is 
uniformly zero throughout the shell. The error introduced is usually negligible, but for 
certain boundary conditions or for sharply varying surface loads, this term may be sig- 
nificant and is consequently retained in this analysis. If 7p ^ and 7p ^ ^ are neglected, 
the problem may be reinterpreted as that of a prestressed but undeformed shell of revo- 
lution perturbed about the undeformed state. With this term retained, the equilibrium 
equation in the meridional direction (eq. (9)) reduces to 
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( 10 ) 


( pN 4) j£ + N {M'5f N e + R^ Q ? + i 

' V? " !( n 4 + N ek« " ? 4 N 4 = ph ‘' 6 

By the same procedure, the remaining equilibrium equations (eqs. (2b) to (2e)) are 
linearized. Solving equations (2d) and (2e) for and and substituting for them 

into equations (2a), (2b), and (2c), eliminates these quantities from the system. The 
parameters defining the geometry of the middle surface can be nondimensionalized by 
using a reference length a, as follows: 



= 1 


r = £ 
a 


(id 


and nondimen sional curvatures can be defined as 




-\ 


k e = — 

R e 


(12) 


Upon completion of these manipulations, the following equilibrium equations result: 

a (S N 4 + rN l,x + n 4«,« - ^ N e) + ^ i m 4 + rk x M |,x ‘ k x | M e + f ( 2k x ‘ k 9 ) M s«,« 

- a|rk x N £ (p £ + 0 + rk x <o^N^J = ra 2 h yU (13a) 


( 2 1 n 4« + rN 49,x + n m) + i|(v * *»)£ - tr 

{'Wf + ls( S 4 + + l( S 4 + s e) V ‘ rk e^ N |e 


M to + |( 3k e “ k x) M ^,x T ~er x e,e 


+ k„M, 


O •• 

= ra^hi/V 


(13b) 
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-"(**"{ + Ve) + I? M « + 2 S M «,x + rM ?,xx 


dr 


_ - — M fl _ + a + 2M fcfl „ + -M fl flfl 

^2 9 dx 9,x r dx f0,0 E,6,x9 r 9,90 

f dr - - dN, 

+ t ‘ rN ^?,x - r it *r *9*9,9 


- & '"4 N I - r l^V r '"s N 4 ,x - ’’{ N «e,e) = ra2h1 '* 


(13c) 


Similarly, the boundary conditions (eqs. (3)) are given as 


= 0 or U = 0 


N ie + ^( 3k <''M M l9 + K N S + N <’)‘ ? '° or v = 0 


M 


l,x rdxl I e ) r 


2M . 


- ^ = 0 or W = 0 


cp^ = 0 or = 0 j 


(14) 


The modified stress-resultant— strain relationships, if physical linearity is assumed, are 


N ? = B^ + fie 0 ) 

N 0 = B ( e 0 + tie i) 

N 40 = B ^ " 


= D ( K £ + IJLK e) 


M, 


= D^/Cg + + D^l - M 2 )«0 


= ‘ ^*$0 


(15) 


J 
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where 


B = 


Eh 

1 - ^ 


D = 


Eh" 


12(l - /i 2 ) 


( 16 ) 


The linearized strain -displacement relationships, from reference 9, reduce to 


s 5 = l( u ,x + V * - *{ W ,x) 




|0"iV2r +_ 2~'ta2 F T ^"^ V ‘ W",ej 


ll'-' 0 v x dr V Vp Vp 

+ -j2. - S£ -f- + — S. k 0 V - — * W 


'r?(v,^ c -w,J 

K -l/dr!ku + ^V W ,00 dr W ,x 
r U + r V r2 dx r 


>(17) 


= a 2 


( 3k x - *fl) n , ( 3k e - k x) y , (^x ~ 3 ^g) dr y W ,xa,dr W , 0 " 
4r 4 > x 4r dx r dx r 2 


The middle -surface rotations are given in terms of displacements as 


n - a( k x U - w ,x) 

= ±(v + S LV.l 

2a \ > x dx r r 




The prestress terms are given nondimensionally as 


(18) 


N| = Be x (x) 


N 0 = B * 0 ( x > 


(19) 
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Reduction to Ordinary Differential Equations 

With equations (13) and (15) to (19), the equilibrium equations can be reduced to 
three partial differential equations with the displacements as the unknown dependent 
variables where the highest order derivative in x is a fourth-order derivative. How- 
ever, since the solution, in general, can only be achieved by numerical techniques, the 
procedure of reference 1 is followed, where dependence on 0 is removed by assuming 
a solution of the separable type and introducing as an additional unknown. This 

procedure yields a set of four second-order ordinary differential equations with variable 
coefficients. The fourth equation is simply the equation for in terms of the dis- 

placements. This reduction in order is essential for the numerical treatment that 
follows. 


A solution is assumed of the form 

U = u(x)(cos n0)e iwt 

V = v(x)(sin n0)e iwt 

W = w(x)(cos n0)e^ a ’^' 

3 

Mj. = m (x)(cos n0)e iCi;t 
^ a 2 x 


( 20 ) 


Defining the perturbation displacements in this manner assures compatibility in the 
0-coordinate. The special case of axisymmetric torsional vibration or torsional 
buckling is precluded in this investigation by the introduction of this form of 0-variation. 
This vibration mode uncouples from the extensional and bending modes (see ref. 14) and 
may be treated directly by interchanging the sines and cosines in equations (20). 


For the buckling problem, the time dependence of the perturbed displacements is 
removed by allowing to (eqs. (20)) to vanish. The perturbed displacements then repre- 
sent a possible equilibrium state. Any loading that would maintain this equilibrium state 
as well as the prestressed equilibrium state is a critical loading for buckling. 


Performing the operations indicated and utilizing the following geometric 
relationships 


dk 

dx r dx V x 


6 = ± ^(K. - k 0 ) 


Mr = - rk x k 0 

dx 2 x 6 


( 21 ) 


j 


14 


which are the Codazzi and Gauss equations, respectively, yields the governing equations, 
as follows: 

F ll u " + G ll u ' + H ll u + G 12 v * + H 12 v + F 13 w " + G 13 w * + H 13 w + G 14 m x’ + H 14 m x = 0 

(22a) 

G 21 u ' + H 21 u + F 22 v " + G 22 v ' + H 22 v + F 23 w " + G 23 w ’ + H 23 w + H 24 m x = 0 (22b) 

F 31 u ” + G 31 u ' + H 31 u + F 32 v " + G 32 v * + H 32 v 

+ F 33 w” + G 33 w’ + H 33 w + F 34 m x " + G 34 m x ’ + H 34 m x = 0 (22c) 

G 41 u' + H 41 u + H 42 v + F 43 w" + G 43 w' + H 43 w + H 44 m x = 0 (22d) 


The same procedure yields the boundary conditions, as follows: 


e H u ' + f j^u + f jgV + e i3 w ' + ^13 W = 0 

or 

u = 

0 

(23a) 

f 21 u + e 22 v ' + f 22 v + e 23 w ' + f 23 w = 0 

or 

v = 

0 

(23b) 

e 31 u ' + f 31 u + e 32 v ' + f 32 v + e 33 w ' + f 33 w + e 34 m x' + f 34 m x = 0 

or 

w = 

0 

(23c) 

f 41 u + e 43 w' = 0 

or 

m x 

= 0 

(23d) 


Primes denote total differentiation with respect to the nondimensional variable x, and 
the coefficients are subscripted for convenience in subsequent matrix manipulations. The 
coefficients F.., G.. , H.. are given in appendix A in terms of the parameters y 


> (24) 


and X, where 


= 1 dr 
y r dx 


X = b 
a 
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o 

For the vibration problem, the frequency parameter O occurs in H 22 , 

and Hgg, where 

= a 2 o) 2 ^(l - u 2 ) 

E 


( 25 ) 


CLOSED-FORM SOLUTION FOR CYLINDER VIBRATIONS 


The vibratory characteristics of a "freely supported" cylinder (simply supported 
but unrestrained in the axial direction) with prestress deformations neglected are well 
known. (See refs. 10, 11, and 12.) Thus, these known results can be used as a check of 
the validity of the governing equations and of the accuracy of the numerical techniques to 
be suggested subsequently in this report. 

When prestress deformations are neglected and the in-plane stresses are constant, 
the equations (22) reduce to ordinary differential equations with constant coefficients 
which, for freely supported boundary conditions, have a solution of the form 


m7TX 


u(x) = A m cos g 


v(x) = B m cos 


w(x) = C m sin 


(26) 


mi tx 


m x (x) = D m sin g 


m = 1,2, . . . 


J 


where S is the length -radius ratio of the cylinder (a = cylinder radius). The classical 
procedure of neglecting prestress deformations to ensure constant coefficients in the field 
equations implies that the cylinder is initially prestressed as a shell with free edges and 
then subsequently supported for vibration. 

Equations (26) are substituted into equations (22) to yield a set of linear homogene- 
ous algebraic equations. For a nontrivial solution to exist, the determinant of the coef- 
ficient matrix of the resultant set of equations must vanish. This procedure leads to the 
characteristic equation 


A 3 (o 6 ) + A 2 (fi 4 ) + A 1 (f2 2 ) + A 0 = 0 
where the coefficients A^ are given in appendix B. 


(27) 
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Equation (27) has been solved for the frequency parameter for an unstressed cir- 
cular cylindrical shell with p = 0.3, A = 0.001, and S = 3; and the results for the axial 
mode m = 1 are given in figure 3. For the purpose of comparing the results of equa- 
tion (27) with those given in references 10, 11, and 12, the in-plane inertia terms are 
dropped from equations (22a), (22b), and (22c), so that equation (27) is reduced from a 
cubic to a linear equation in In all studies performed, the frequencies calculated, 

when in-plane inertias were neglected, agreed closely with those given in references 10, 
11, and 12, and all trends observed in these references were similarly verified. 

The effect on the natural frequencies due to the neglect of in-plane inertia is illus- 
trated in figure 3. The error introduced by this approximation decreases as n 
increases, the largest error occurring at n = 2. For n larger than 5 the error is 
negligible. For n = 0, the fundamental frequency, which corresponds to pure torsional 
oscillations, is excluded when the assumption of negligible in-plane inertias is made. 



Figure 3.- Natural frequencies of freely supported circular cylindrical shell. 

X = — = 0.001; S = 4 = 3; m = 1. 
a a 
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DEVELOPMENT OF NUMERICAL PROCEDURE 


Development of Difference Equations 

A numerical procedure is needed for those shells of revolution and loading condi- 
tions which do not admit a solution in closed form. The meridian of the shell is divided 
into increments and a three -point central difference method is used to reduce the differ- 
ential equations to algebraic form. The distance measured along the meridian between 
adjacent stations is constant and is represented nondimensionally by A where 

i = = | < 28 > 

and where the subscript i on symbols and matrices indicates the evaluation of the sub- 
scripted variable or matrix at the ith station, i = 0, 1, 2, . . ., N and where 

S total meridional arc length of the nondimensional shell, s/a 

N total number of intervals 


The three -point difference formulas, when applied at the ith station for some func- 
tion z(x), are 


z i" * T2 ( z i-l “ 2z i + z i+l) 


Z i' *2^-1 + Z i + l) 


(29) 


Reference 2 indicates that this simple approximation leads to sufficiently accurate results. 
The governing equations (22) may be written in matrix form at station i as 


F i Z i M + G i Z i’ + H i Z i = 0 

where 


F 11 

0 

F 13 

0 

0 

F 22 

F 23 

0 

F 31 

F 32 

F 33 

F 34 

0 

0 

F 43 

0 


i 


(30) 


(31a) 
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and 


G 11 

G 12 

G 13 

G 14 

G 21 

G 22 

G 23 

0 

G 31 

G 32 

G 33 

G 34 

G 41 

0 

G 43 

0 

*11 

H 12 

H 13 

H 14 _ 

H 21 

H 22 

H 23 

H 24 

H 31 

H 32 

H 33 

H 34 

h 4 i 

H 42 

H 43 

H 44_ 


fu 1 





in 


Similarly, all the boundary equations (23) may be written in matrix form 
boundaries, i = 0 and i = N, as 


where 


a o e o z o' + (“o f o + ^o) z o = ° 

a N e N Z N + (“N f N + %) Z N = 0 


e ll 

0 

e 13 

0 


0 

e 22 

e 23 

0 


e 31 

e 32 

e 33 

e 34 


0 

0 

e 43 

0 

0,N 


(31b) 


(31c) 


(3 Id) 


at the 


(32) 


(33a) 
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(33b) 


and where 


[ 0,N 


“0,N 


^0,N = 


f ll 

f 12 

f 13 

0 

f 21 

f 22 

f 23 

0 

f 31 

f 32 

f 33 

f 34 

£41 

0 

0 

0 

‘'ll 

0 

0 

0 

1 

a 22 

0 

0 

i 

0 

<*33 

0 

1 

0 

0 


'11) 

0 

0 

0 

e 

" “22) 

0 

0 


0 

( l - “33) 

0 


0,N 


! o,n 


( x - “44) 


(34) 


0,N 


The elements 


take on the value 1 or 0 depending on the prescribed conditions. 


The a- and /3-matrices (eqs. (34)) are used to select the prescribed boundary con- 
ditions. If, for example, u=0 is prescribed at i=0 then =0 and if u is 

not prescribed at i = N (i.e., if the u displacement is unrestrained in the meridional 
direction of the undeformed shell), then ^ a ll) = L If desired, the present theory can 

be extended to allow for elastic and directional supports in the boundary conditions by 
appropriate redefinition of the a- and /3-matrices. 


When equations (29) are applied to equations (30) and (32), the governing equations 
become 




z. 


i-1 


+ 




+ 




= 0 


(i = 0, 1, 2, . . ., N - 1, N) (35) 


20 



and the boundary equations become 


«0 e 0 / \ a n e o 

(“o f o + 0o) z o + -^ z i = ° 

2i~ Z N-1 + (“N*N + %) Z N + 2 A N Z N+1 




a N e N 


Z-kt. 1 — 0 


(36) 


Equations (36) can be solved for Z_^ and Z N+ j and the results can be substi- 
tuted into equation (35) to yield, at i = 0, 


,-1 


2F, 


fOp) 2° | IH il| + f 

2A a 2 ' 2A 0 " A 2 +f 0 


0) 


and, at i = N, 


+ ftu z n + 


“0 e 0 


oro T 2^ 


A 2 2A a 2 2A 


N 


!n/2n + M V h 

2AU2 2a) % 


2F- 




+ %( Z N 


Zj = 0 

(37) 


where 


a N e N 


2A 


Zn + % 

, a2 2A 


-1, 


_N 

,A2 


_N ' 
2A i 


+ I 


Z N-1 " 0 


(38) 


I = 


0 

1 

0 

0 


0 

0 

1 

0 


The difference equations (35), (37), and (38) constitute a complete set of field equations 
governing the behavior of the perturbed state. 


Numerical Solution 

The problem is now one of solving a set of homogeneous equations (eqs. (35), (37), 

and (38)). This set constitutes an eigenvalue problem such that the mode shape Z i is 

9 1 

the eigenvector and the frequency parameter SI* is the corresponding eigenvalue for the 
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vibration problem. The eigenvalue for the buckling problem is contained in all barred 
terms. The fourth equation, being simply the definition of m x in terms of displace- 
ments, will not contain an eigenvalue for either problem. For a nontrivial solution to 
exist, the determinant of the coefficient matrix must vanish. 

The coefficient matrix will be a 12-element -wide band matrix. A convenient tech- 
nique for solution of this problem can be formulated by modifying the method of refer- 
ences 7 and 8. Such a modification is presented herein to handle the free vibrations and 
buckling of a prestressed shell governed by four second-order difference equations with 
two-point boundary conditions. 

Define the following (4 x 4) matrices: 


F. G- 
^ A 2 2A 

2F. 


•'v 


B. = H. - 


A2 


F. G. 
C. = + — 

1 A 2 2A 


D 0 =O! 0 


-1 


e 0/ F 0 « 0 \ (rj 

Sill* “SI H o 


2F o\ 

-J- +f o 


+ £r 


d N=«n 


"!n/^ + £nV 1 | 


2F 


2A\A2 T 2A j p" A 


N +f 

2 N 


% 


En — ' 

0 0 2A 


!o.£oV7fo + S\ + I 

a2 2A a 2 2A 


N 

Ext = -«AT ~ 

N N 2A 


i A 2 2A 


- 1 , 


f n_2n) 

2A I 


+ I 


(39) 


Equations (35), (37), and (38) may now be written as 

+ B^Z^ + C^Zj + ^ = 0 (i = 1, 2, . . ., N - 1) 


D o z o + E 0 Z 1 “ 0 


e n z n-i + d n z n = 0 


(40) 

(41) 

(42) 


L_ 
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For such a set of homogeneous equations, a recursion formula for may be 

written as 

(i = 1, 2, . . ., N - 2) 
and 

(i=0 if Z Q ^0) 

(i = N - 1 if Z N t 0) 


Z i + P i Z i + l " 0 


( 43 ) 


where Pj is a (4 x 4) recursion matrix. To find P^, combine equation (43) and equa- 
tion (40) to obtain 


,-1 


Z i + ( B i - A i P i-l)" C i Z i+l = 0 (i = 1, 2, . . N - 1) (44) 


Comparison of equation (44) with equation (43) shows that 

.-1 


P i = ( B i " A i P i-l)’ C i (i = 1, 2, . . N - 1) (45) 

Comparison of equation (41) with equation (43), the latter written at i = 0, shows that 


p 0 = D o E o 


(46) 


From equations (45) and (46), Pj may be found at all points with the exception of the 
point i = N. This process of determining all required values of P^ in terms of Pq 
is in essence a Gaussian elimination process. 

Equation (43) written at i = N - 1, in combination with equation (42), yields 

( d n- e n p n-i) z n = ° < 47 > 


If Z^j ^ 0, then for a solution to exist, 


°N ' % P N-1 


= 0 


(48) 


Therefore, any frequency parameter (or buckling parameter in the corresponding 

stability problem) which satisfies equation (48) contains a natural frequency (or critical 
load) of the system. The natural frequencies can be found by trial and error by selecting 
successive values for calculating the matrices of equations (39), and using equa- 
tions (45) and (46) to evaluate the determinant in equation (48). This procedure is con- 
tinued until the desired zeroes of the determinant are found. 

The method must be slightly modified for the case Z^ = 0. Substituting equa- 
tion (43) written at N - 2 into equation (40) written at N - 1 yields 
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( 49 ) 


( B N-1 " A N-1 P N-2) Z N-1 _ 0 
If Z N j = 0, then Zj = 0 from equation (43), and the solution is trivial. Therefore 


B N-1 " A N-l P N-2 " 0 


(50) 


Consequently, for the case Z^ = 0, equation (50) is used in place of equation (48) in the 
elimination process. 


After the natural frequencies have been found, the mode shapes are determined by 
first solving equation (47) for Z^ in terms of an appropriate normalizing factor. For 
the case where Z-^ = 0, equation (49) is used to solve for Z^^. The remaining Z^'s 
are then determined by using the recursion formula, equation (43). 


This numerical procedure is particularly well suited for use with a large number 
of stations since only the band elements need be retained during the computation process. 
References 1, 2, and 4 give further advantages in using this general method of solution. 


In an investigation of buckling, the procedure must be slightly altered. The exter- 
nal axisymmetric loading P ^ and P can be given by 


P^(x) = Kp^(x) 
P(x) = Kp(x) 


(51) 


where and p characterize the form of the external applied load and K is a con- 
stant governing its magnitude. If linearity between the external applied loads and the 
prebuckling in-plane stresses and deformations is tacitly assumed, the prebuckling quan- 
tities may be written as 


and it follows that 


"x = KE x 


e o ~ 


cp ^ = K<1>£ 


de. 


= K 


d_E^ 
dx 

dx dx 


x 

dx 

de 


d <p t d<£ t 

_1= K — £ 
dx dx 


J 


(52a) 


(52b) 
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where E x , E^, characterize the in-plane stress and deformation due to 

and p. These quantities can be found, for example, from the stress program of refer- 
ence 1. Thus, K becomes the buckling load parameter, and any value of K which 
allows two (or more) equilibrium states to exist simultaneously defines a critical 
buckling load for the system. The numerical procedure is followed directly, as before, 
with trial values of K rather than selected to satisfy equation (48) or equation (50). 

With K known, equation (51) yields the critical loading for the particular harmonic 
wave number n investigated. 

Although equation (48) (or eq. (50)) contains all the roots of the system of equa- 
tions (40), (41), and (42), this method of elimination introduces spurious singularities in 
the determinants of equation (48) or (50). For some shell configurations and boundary 
conditions, it is found that the roots and singularities very nearly coincide, and the usual 
predictor -corrector methods fail to indicate a root if the increments given to the fre- 
quency parameter or buckling load parameter are too large. Moreover, some of these 
singularities are associated with a change in sign in the value of the determinant even 
though no zero exists at that value of the frequency or buckling parameter. A technique 
for avoiding this difficulty is presented in reference 15. 

EXAMPLES OF APPLICATION OF NUMERICAL PROCEDURE 


Cylinder Vibrations and Buckling 


An indication of the accuracy of the numerical procedure is obtained by investi- 
gating the free vibration and buckling of an unstressed, freely supported, cylindrical 
shell. The freely supported boundary conditions are introduced by defining 


10 0 0 


a 


0,N 


0 0 
0 0 


0 

0 


0 0 0 


0 


The particular shell investigated has the parameters p = 0.3, X = 0.001, and S = 3, 
where the reference length is taken to be the cylinder radius, so that X is the ratio of 
the thickness to the radius and S is the ratio of the length to the radius. The calcula- 
tions are based on 200 intervals. 

The lowest frequency for each circumferential harmonic mode number is found for 
an unstressed shell and compared with the corresponding results calculated using the 
closed-form solution of equation (27) with m taken as 1. If S is held constant and 
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m is varied for any n, it can be shown that the lowest frequency will result for m = 1. 
The n um erical and exact results agree to three significant figures or better for all values 
of n from 1 to 15. 

The same cylinder is investigated with one edge clamped and one edge freely sup- 


ported. The clamped edge is introduced by defining 



o 

0 

0 

0 


0 

0 

0 

0 

“0 = 

0 

0 

0 

0 


0 

0 

0 

1 


A comparison between the two edge conditions for the lowest frequency for each circum- 
ferential harmonic mode number is presented in figure 4. As expected, introducing the 

clamped edge increases the nat- 
ural frequencies of the system, 
the lowest natural frequency 
having a relative increase based 
on the freely supported cylinder 
of 26 percent. As n increases, 
the effect of the edge conditions 
appears to diminish rapidly. 

The buckling load is found 
for a closed circular cylinder 
subjected to uniform hydrostatic 
pressure. The results from the 
numerical procedure are com- 
pared with the closed -form solu- 
tions based on Donnell theory 
given in reference 13. Pre- 
buckling deformations are neg- 
lected in the derivation of refer- 
ence 13; consequently, to allow 
for consistent comparisons, the 
assumption is made that 7p ^ = 0. 

For a closed cylinder under 

- ^9 

external pressure, = — . 

Thus, unit stress conditions 

a 



Circumferential harmonic wave number, n 
Figure 4.- Comparison between the natural frequencies of a circular cylinder 
freely supported at both ends and freely supported at one end and clamped 

at the other end. A = 0.001: S = - = 3: u = 0.3. 
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characterizing this loading are introduced by letting Eg = -1 and E x = - j. Hence, 
from equations (19) and (52) 



where Ng may be taken to be Pp, P is the uniform external pressure (P is a nega- 
tive quantity), and p is the radius of the cylinder measured to the middle surface. 

Thus, the buckling load parameter is 


-(l - h 2 )p 
EX 


(54) 


The corresponding K from reference 13 in terms of the nondimensional param- 
eters used in this development is given by 


A 2 

” 1 + (t) 2 ] 

2 

s(i- 

m 2 ) 

12S 2 

[Mtf 

1 

Miff 

iHff 


(55) 


The lowest buckling load occurs for this shell at n = 9. The closed-form solution is 
generally accepted as giving accurate results for n greater than 4; thus, the closed-form 
value for K = 0.888 x lO - ^ at n = 9 is a valid basis for comparison. By use of the 
numerical procedures outlined herein, the critical buckling load was similarly found to 
be at n = 9 and gave a relative error of 0.68 percent. For n = 1 and n = 2 it is 
known that the lowest buckling load calculated, when the closed-form solution is used, 
will be in error on the high side. The numerical procedure gave relative errors based 
on the closed-form solution of -18.13 percent for n = 1 and -14.23 percent for n = 2. 
These errors indicate that for these cases the use of the numerical procedure coupled 
with the more consistent shell theory yields more accurate results. 


Vibrations of a Shell of Positive Gaussian Curvature 

The previous examples are applications of the numerical procedure for both 
buckling and vibration of systems governed by equations with constant coefficients. The 
capabilities of the numerical procedure are further demonstrated by investigating the 
vibration characteristics of a shell of positive Gaussian curvature with a constant positive 
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meridional curvature and comparing them with those of a corresponding shell of zero 
Gaussian curvature (cylinder). The correspondence is introduced by requiring the total 
surface area (and thus the mass) of both shells to be equal. Furthermore, the meridional 
lengths of both shells are selected to be equal and thus the perpendicular distance from 
the axis of revolution to the centroid of the meridional curve is the same for both shells. 
This distance is selected as the reference length for each shell. 

The shapes of both shells investigated are illustrated in figure 5 along with a tabu- 
lation of the equations used in the numerical program to define the geometries. The 
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Figure 5.- Geometry of shells investigated. 



lowest natural frequencies for successive circumferential harmonic mode numbers, 
based on freely supported edge conditions, have been found for each unstressed shell. 
The results are presented in figure 6. 



Circumferential harmonic wave number, n 
Figure 6.- Comparison between frequencies of shell of positive Gaussian 
curvature and a corresponding cylinder 


The lowest frequency for the positive -curvature shell was considerably higher than 
the corresponding frequency for the zero curvature shell for each circumferential mode 
number greater than 2. The maximum increase tabulated was in excess of 1,000 percent 
indicating a large increase in stiffness due to the curvature of the meridian. 

CONCLUDING REMARKS 

A set of linear equations governing the infinitesimal vibrations and buckling of 
axisymmetrically prestressed shells has been developed and both the in-plane inertia 
and prestress deformation effects have been retained in the analysis. The equations 
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derived are consistent with first-order thin-shell theory and can be used to describe the 
behavior of shells with arbitrary meridional configuration having moderately small pre- 
stress rotation. 

A numerical procedure has been given for solving the governing equations for the 
natural frequencies, buckling loads, and associated mode shapes for a general shell of 
revolution with homogeneous boundary conditions. The numerical procedure uses matrix 
methods in finite -difference form coupled with a Gaussian elimination to solve the gov- 
erning eigenvalue problem. 

Examples of applications of the numerical procedure have been presented. Results 
for natural vibration frequencies and buckling under hydrostatic pressure of a simply 
supported cylinder were found to be in good agreement with previously published results. 
A brief comparison between the vibration characteristics of two geometrically similar 
shells having positive and zero Gaussian curvature displayed the markedly greater stiff- 
ness and resulting higher frequencies of the doubly curved shell. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., September 14, 1966, 

124-08-06-11-23. 
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APPENDIX A 


COEFFICIENTS OF EQUATIONS (22) AND (23) 

Coefficients of Equations (22) 

The coefficients of terms in the governing equations (22) are defined as follows: 


F 11 = 1 


F 13 " 


F 22 “ + ~ ^96 " k x ) 2 + |( e x + e e) 

F 23 = x2<1 2 :/ )n ( 3k 8 ■ “*) 


F 31 “ F 13 


F 32 " F 23 


33 


F 34 = V 


X 2 (1 - n) 

2n 2 

12 

[r 2 


■ + (l + ii)r z 


+ e x+ cp L 




F 43 ~ F 34 


G U-r 


G i2 = + - k e) ( 3k e - k x) - £(*x + *e) 




G 14 = X 2 (l -M 2 )k x 


G 21 = _G 12 




_ _ _ 2 
+ k x e x - ~ MV* + Mi 


G 22‘ il f Jil r-^ JiX ( 3k e- k ^ 


2 ^ - >'( 5k * - 3k e) 


4*r( e x +e e) + 4( e *'" e 8) 
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APPENDIX A 


r - a (1 - M) n 
G 23 “ 24r 


dk 

2(l +M )rk e --^+3y(k x -k 0 ) 


+ {l+jjjn - 

2r ? 


- „ d< ?f 

G 31 = _G 13 " 2y< ^ " 2 “5“ 


G . X (1 - M)n 
32 2 4r 


dkx 

Syk x - yk e (5 + 2 ja) - — 


dx 


- ^-(1 + 

2r ? 


, _ X Z (1 - p.) 

'33 - - 12 


(l + M)(2yk x k 0 + r 3 ) + 


de x _ 2 _ d <P, 

+ ye x + -5T +y<p ^ +2 ^1I 


G 34 = x 2 ( 2 - M)(l - P 2 )y 


G 41 - -G 14 


G 43 = a2 (i " M 2 )^y 

H U = -(***<, - '>' 2 - {L ^~ - i % Jli 


' S2 ( 5 * + *«) + 


dk~ 

(1 - v)ykx + _ 


8r^ 

,2-2 , d fl 0 2 
'V k x * { +k x-5f + 


V^x 


H 


12 " “ 2r 


X 2 (l - n)n y 

■(3kx-k 0 )( 3k e -k x ) 

12r 

8 


+ (1 + jtxJkjjkg 


- g(s + •») 




dk. 


13 
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X 2 (l - n)n 2 y 

( 3k x - k e) 

12r 2 
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— (1 - u)n 2 — 
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h 14 = x 2 (i - /x 2 )(i - fib** 
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' - x2(1 12 7 )n [» + chv, * £(<*,*9 - 7k x 2 - 3k e 2 ) - i^( 5k e - 3k x) 
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(1 + fi)n 2 

12 

r 2 


k k. 


+ 1 s( e x + e e) - k / e e - jVi( , i + e e) + iLz ^ bL { k x * 


- ko 2^2 + n 2 


v'i 




® dx 


dk. 


r -gf - 2 - 2(1 + j* )n 


k e * ( 3k e - 1 S.)(’' 2 + k x k e) 


- k^e. + -/ >1>r V, - *H^kJp 2 + 

r 0 0 2r ? 2r y s 2r dx 


H 2 4 = . kfl 


«3i = -r( k e + " k x) + - ^r(y\ - r ^ * *\\) 


n“ 

2r 2 


dk^ 

yk * - ^ - 3 IT; 


dk. 




de. 


(l-M)n 2 k 2 

2r2 


<Pi 


-(’' k x + ^)'‘ > { 2 -( 2 k x' , ’{ +f ‘ y )S 


h 32 “ - ?( k e + “ k x) + x2(1 24/ )n 2(1 - - Ac + 2 A - ^5 k o) + y -£T 


+ 3y ( k e- k x) + k x k e( 3k 9- k x) 


■ n k S + (1 - M)ny ? q - Jib k ? 2 . ^ gj 

r & & 2r s 2r ® s r dx 


33 


I 


APPENDIX A 


h 33 - V - We - k e 2 - - <t 12 'i >n2 
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Coefficients of Equations (23) 

The coefficients of the terms associated with the boundary conditions (eqs. (23)) are 
defined in the following equations: 

e H= 1 

e 13 = "^5 


e 22 " ^f^+ X (1 96 ' ^(3k 0 - k x ) 2 + ±(e x + e Q ) 
e 23 - - k x) 


e 31 “ e 13 
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COEFFICIENTS OF EQUATION (27) 
The coefficients of equation (27) are defined as follows: 


A 0 = 


a ll a 22 a 33 " a l 1 ( a 23 ) 2 " a 22( a 13) 2 ' a 33( a 12 ) 2 + 2a 12 a 23 a 13 


*44 
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